setwd("~/Nahdha Replication MedPol/Figures")




####Figure 1####
secular.mean<-c( 2.731,2.629, 2.846,  3.033, 3.236)
secular.label<-c( 'Strongly \n Disagree','Somewhat \n Disagree','Neither Agree \n nor Disagree','Somewhat \n Agree','Strongly \n Agree' )
par(mar=c(5, 4, 4, 4))
barplot(secular.mean, main="Secularism",
        ylab="Mean Favorability",names.arg = secular.label, ylim=c(2,4),xpd = FALSE,col='gray25')
education.mean<-c( 2.525 , 2.933,  3.015, 3.234, 3.777)
education.label<-c( 'No formal \n education','Elementary\n','High School\n','University\n','Higher\n Education' )


barplot(education.mean, main="Education",
        ylab="Mean Favorability",names.arg = education.label, ylim=c(2,4),xpd = FALSE,col='gray75')

####Figure 2####
library(ggplot2)
library(stringr)
data.education<- read.csv("q3.csv")
head(data.education)


model1Frame <- data.frame(Education = data.education$education,
                          Coefficient = data.education$difference,
                          SE = data.education$se,
                          modelName = "Education") 

model1Frame$Variable <- factor(model1Frame$Education, as.character(model1Frame$Education))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(model1Frame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2), color='black')
zp1 <- zp1  + theme_bw(base_size = 14)+ theme(legend.position="none") + xlab("Highest Level of Education Attained") + ylab("Change in Favorability") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().

data.rural<- read.csv("rural.csv", sep="\t")
head(data.rural)

data.rural<-data.rural[-c(),]
model1Frame <- data.frame(Education = data.rural$rural,
                          Coefficient = data.rural$difference,
                          SE = data.rural$se,
                          modelName = "Rural") 

model1Frame$Variable <- factor(model1Frame$Education, as.character(model1Frame$Education))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(model1Frame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2), color='black')
zp1 <- zp1  + theme_bw(base_size = 14)+ theme(legend.position="none") + xlab("") + ylab("Change in Favorability") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


secular<- read.csv("q8.csv")
head(secular)


model1Frame <- data.frame(Education = secular$secular,
                          Coefficient = secular$difference,
                          SE = secular$se,
                          modelName = "Secular") 

model1Frame$Variable <- factor(model1Frame$Education, as.character(model1Frame$Education))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(model1Frame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2), color='black')
zp1 <- zp1  + theme_bw(base_size = 14)+ theme(legend.position="none") + xlab("Religion and Politics should be Separate") + ylab("Change in Favorability") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().


service<- read.csv("q13.csv")
head(service)


model1Frame <- data.frame(Education = service$service,
                          Coefficient = service$difference,
                          SE = service$se,
                          modelName = "Secular") 

model1Frame$Variable <- factor(model1Frame$Education, as.character(model1Frame$Education))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(model1Frame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2), color='black')
zp1 <- zp1  + theme_bw(base_size = 14)+ theme(legend.position="none") + xlab("Public Service has Improved since 2011") + ylab("Change in Favorability") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().

party<- read.csv("q18.csv", sep="\t")
head(party)
party<-party[-c(5),]

model1Frame <- data.frame(Education = party$party,
                          Coefficient = party$difference,
                          SE = party$se,
                          modelName = "Party") 

model1Frame$Variable <- factor(model1Frame$Education, as.character(model1Frame$Education))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(model1Frame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2), color='black')
zp1 <- zp1  + theme_bw(base_size = 14)+ theme(legend.position="none") + xlab("Which political party do you support most?") + ylab("Change in Favorability") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 
print(zp1) # The trick to these is position_dodge().

####Figure 3 Plots####



conf=.95
title="No Formal Education"
xlabel=""
ylabel="Change in Favorability"
factor_labels=c("Baseline","Reform Prime")


# Get coefficients of variables
beta_1 =  .7469616  
beta_3 =   -.6757209  

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =     .03185781  + (x_2^2)*   .05913456   + 2*x_2*   -.0312186

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)


conf=.95
title="Elementary Education"
xlabel=""
ylabel="Change in Favorability"
factor_labels=c("Baseline","Reform Prime")


# Get coefficients of variables
beta_1 =  .2588769  
beta_3 =  -.2718607  

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =    .01288147  + (x_2^2)*    .02719556  + 2*x_2*   -.01249245

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)




conf=.95
title="Secularism: Neither Agree nor Disagree"
xlabel=""
ylabel="Change in Favorability"
factor_labels=c("Baseline","Reform Prime")


# Get coefficients of variables
beta_1 =  .4428872  
beta_3 =   -.6588398 

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =    .0241985 + (x_2^2)*     .06092094   + 2*x_2*    -.02482854

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)



